#Load the data
data(UCBadmit)
d <- UCBadmit
d
## dept applicant.gender admit reject applications
## 1 A male 512 313 825
## 2 A female 89 19 108
## 3 B male 353 207 560
## 4 B female 17 8 25
## 5 C male 120 205 325
## 6 C female 202 391 593
## 7 D male 138 279 417
## 8 D female 131 244 375
## 9 E male 53 138 191
## 10 E female 94 299 393
## 11 F male 22 351 373
## 12 F female 24 317 341
summary(d)
## dept applicant.gender admit reject applications
## A:2 female:6 Min. : 17.00 Min. : 8.0 Min. : 25.0
## B:2 male :6 1st Qu.: 45.75 1st Qu.:188.2 1st Qu.:291.5
## C:2 Median :107.00 Median :261.5 Median :374.0
## D:2 Mean :146.25 Mean :230.9 Mean :377.2
## E:2 3rd Qu.:154.00 3rd Qu.:314.0 3rd Qu.:452.8
## F:2 Max. :512.00 Max. :391.0 Max. :825.0
#Fit the model with map2stan
#First, creat a dummy variable for "application.gender"
d$male <- ifelse( d$applicant.gender=="male" , 1 , 0 )
#Second, make index for each dpartment
d$dept_id <- coerce_index( d$dept )
#Third, get rid of "." in the colum name "application.gender"
colnames(d) <- sub(".","_",colnames(d),fixed = TRUE)
mQ1.stan <- map2stan(
alist(
admit ~ dbinom( applications , p ) ,
logit(p) <- a_dept[dept_id] + bm*male ,
a_dept[dept_id] ~ dnorm(0,10) ,
bm ~ dnorm(0,10)
) , data=d, chains = 4)
## In file included from file10a4d398053.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/var.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/tools/config.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config.hpp:39:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config/compiler/clang.hpp:196:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command line>:6:9: note: previous definition is here
## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1
## ^
## 1 warning generated.
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.050203 seconds (Warm-up)
## 0.038984 seconds (Sampling)
## 0.089187 seconds (Total)
##
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.048116 seconds (Warm-up)
## 0.044374 seconds (Sampling)
## 0.09249 seconds (Total)
##
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.051587 seconds (Warm-up)
## 0.05759 seconds (Sampling)
## 0.109177 seconds (Total)
##
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 4).
##
## Chain 4, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.056264 seconds (Warm-up)
## 0.04476 seconds (Sampling)
## 0.101024 seconds (Total)
##
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 3.8e-05 seconds (Sampling)
## 4.1e-05 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
#test the model
precis(mQ1.stan,depth = 2)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_dept[1] 0.68 0.10 0.51 0.83 2799 1
## a_dept[2] 0.64 0.12 0.46 0.83 3026 1
## a_dept[3] -0.58 0.07 -0.69 -0.46 4000 1
## a_dept[4] -0.61 0.09 -0.74 -0.47 4000 1
## a_dept[5] -1.06 0.10 -1.22 -0.91 4000 1
## a_dept[6] -2.64 0.16 -2.89 -2.38 4000 1
## bm -0.10 0.08 -0.23 0.03 2322 1
pairs(mQ1.stan)
summary(mQ1.stan)
## Inference for Stan model: admit ~ dbinom(applications, p).
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75%
## a_dept[1] 0.68 0.00 0.10 0.49 0.62 0.68 0.75
## a_dept[2] 0.64 0.00 0.12 0.41 0.56 0.64 0.72
## a_dept[3] -0.58 0.00 0.07 -0.73 -0.63 -0.58 -0.53
## a_dept[4] -0.61 0.00 0.09 -0.79 -0.67 -0.61 -0.56
## a_dept[5] -1.06 0.00 0.10 -1.26 -1.13 -1.06 -0.99
## a_dept[6] -2.64 0.00 0.16 -2.96 -2.74 -2.63 -2.52
## bm -0.10 0.00 0.08 -0.26 -0.16 -0.10 -0.05
## dev 96.13 0.09 3.80 90.78 93.25 95.48 98.26
## lp__ -2597.29 0.04 1.90 -2601.92 -2598.35 -2596.96 -2595.85
## 97.5% n_eff Rhat
## a_dept[1] 0.89 2799 1
## a_dept[2] 0.87 3026 1
## a_dept[3] -0.44 4000 1
## a_dept[4] -0.44 4000 1
## a_dept[5] -0.87 4000 1
## a_dept[6] -2.34 4000 1
## bm 0.05 2322 1
## dev 105.38 1917 1
## lp__ -2594.61 1917 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jan 5 18:19:29 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
#Fit the model by brms, I'm confused by how to create the outcome variable, so I tried different ways
#(in chapter10, 10.1.4, an aggregated binomial uses cbind to build the outcome variable)
mQ1.brms1 <- brm(cbind(admit,reject) ~ 0 + dept + applicant_gender,#separate intercept for each dept
#family = "binomial",#Multivariate models are not yet implemented for family 'binomial'.
prior=set_prior("normal(0,10)", class="b"),
data = d)
## Compiling the C++ model
## Start sampling
summary(mQ1.brms1)
## Family: gaussian (identity)
## Formula: cbind(admit, reject) ~ 0 + dept + applicant_gender
## Data: d (Number of observations: 12)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## admit_deptA 2.11 10.04 -17.51 21.77
## admit_deptB 1.32 10.06 -18.32 21.73
## admit_deptC 0.04 9.86 -18.69 19.58
## admit_deptD 0.08 10.10 -19.49 20.17
## admit_deptE -0.38 10.09 -19.65 19.14
## admit_deptF -1.38 9.62 -20.19 17.57
## admit_applicant_gendermale 2.14 9.96 -17.45 21.62
## reject_deptA -0.43 9.73 -19.59 18.84
## reject_deptB -0.26 9.65 -19.48 18.51
## reject_deptC 0.99 10.05 -18.69 20.78
## reject_deptD 0.89 9.77 -19.19 20.29
## reject_deptE 1.03 10.02 -18.70 21.04
## reject_deptF 1.89 9.83 -16.94 21.22
## reject_applicant_gendermale 1.47 10.12 -18.04 21.42
## Eff.Sample Rhat
## admit_deptA 4000 1
## admit_deptB 4000 1
## admit_deptC 4000 1
## admit_deptD 4000 1
## admit_deptE 4000 1
## admit_deptF 4000 1
## admit_applicant_gendermale 4000 1
## reject_deptA 4000 1
## reject_deptB 4000 1
## reject_deptC 4000 1
## reject_deptD 4000 1
## reject_deptE 4000 1
## reject_deptF 4000 1
## reject_applicant_gendermale 4000 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma(admit) 203.77 39.88 143.53 298.03 4000 1
## sigma(reject) 256.81 52.83 178.56 377.04 4000 1
## rescor(admit,reject) 0.61 0.17 0.19 0.86 4000 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(mQ1.brms1, ask = FALSE)
#then, I create a outcome variable for the model, calculate the probablity of admit in each department
d$admit_p <- d$admit/d$applications
mQ1.brms2 <- brm(admit_p ~ 0 + dept + applicant_gender,#separate intercept for each dept
prior=set_prior("normal(0,10)", class="b"),
data = d)
## Compiling the C++ model
## Start sampling
summary(mQ1.brms2)
## Family: gaussian (identity)
## Formula: admit_p ~ 0 + dept + applicant_gender
## Data: d (Number of observations: 12)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## deptA 0.74 0.09 0.60 0.89 4000 1
## deptB 0.67 0.08 0.53 0.82 3748 1
## deptC 0.37 0.08 0.23 0.52 3762 1
## deptD 0.36 0.08 0.20 0.50 3424 1
## deptE 0.28 0.08 0.13 0.43 3135 1
## deptF 0.08 0.08 -0.07 0.24 2632 1
## applicant_gendermale -0.04 0.06 -0.14 0.08 2942 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 0.09 0.06 0.04 0.21 483 1.01
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(mQ1.brms2, ask = FALSE)
#Test whether the coefficient for gender is different from 0 in the brms model
hypothesis(mQ1.brms1,"admit_applicant_gendermale = 0")
## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (admit_applicant_... = 0 2.14 9.96 -17.45 21.62 NA
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.
plot(hypothesis(mQ1.brms1,"admit_applicant_gendermale = 0"))
hypothesis(mQ1.brms2,"applicant_gendermale = 0")
## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (applicant_gender... = 0 -0.04 0.06 -0.14 0.08 NA
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.
plot(hypothesis(mQ1.brms2,"applicant_gendermale = 0"))
#fit model with interactions between department and gender
#(1)
mQ1.brms3 <- brm(cbind(admit,reject) ~ 0 + dept * applicant_gender,
prior=set_prior("normal(0,10)", class="b"),
data = d)
## Compiling the C++ model
## Start sampling
summary(mQ1.brms3)
## Family: gaussian (identity)
## Formula: cbind(admit, reject) ~ 0 + dept * applicant_gender
## Data: d (Number of observations: 12)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## admit_deptA 2.06 10.15 -17.47 21.87
## admit_deptB 1.23 10.16 -18.30 20.62
## admit_deptC -0.10 10.09 -19.59 19.70
## admit_deptD -0.09 10.03 -19.64 19.18
## admit_deptE -0.25 10.18 -19.82 19.41
## admit_deptF -1.31 9.96 -20.82 18.62
## admit_applicant_gendermale 2.19 9.69 -16.89 20.69
## admit_deptB:applicant_gendermale 1.06 9.94 -18.19 20.13
## admit_deptC:applicant_gendermale -0.06 9.96 -18.90 19.25
## admit_deptD:applicant_gendermale 0.08 9.56 -18.49 18.82
## admit_deptE:applicant_gendermale 0.02 10.16 -19.75 19.80
## admit_deptF:applicant_gendermale -0.70 9.89 -19.97 18.17
## reject_deptA -0.61 9.62 -19.22 17.86
## reject_deptB -0.36 9.88 -19.93 18.91
## reject_deptC 1.07 10.08 -18.49 20.72
## reject_deptD 1.01 9.82 -18.33 20.24
## reject_deptE 0.86 9.93 -18.64 20.32
## reject_deptF 1.79 9.79 -17.66 20.80
## reject_applicant_gendermale 1.62 9.88 -17.36 20.90
## reject_deptB:applicant_gendermale -0.05 9.83 -19.20 19.27
## reject_deptC:applicant_gendermale 0.23 10.08 -19.14 20.00
## reject_deptD:applicant_gendermale 0.54 9.54 -18.35 19.03
## reject_deptE:applicant_gendermale 0.42 9.51 -18.18 19.38
## reject_deptF:applicant_gendermale 1.01 10.00 -18.94 20.76
## Eff.Sample Rhat
## admit_deptA 4000 1
## admit_deptB 4000 1
## admit_deptC 4000 1
## admit_deptD 4000 1
## admit_deptE 4000 1
## admit_deptF 4000 1
## admit_applicant_gendermale 4000 1
## admit_deptB:applicant_gendermale 4000 1
## admit_deptC:applicant_gendermale 4000 1
## admit_deptD:applicant_gendermale 4000 1
## admit_deptE:applicant_gendermale 4000 1
## admit_deptF:applicant_gendermale 4000 1
## reject_deptA 4000 1
## reject_deptB 4000 1
## reject_deptC 4000 1
## reject_deptD 4000 1
## reject_deptE 4000 1
## reject_deptF 4000 1
## reject_applicant_gendermale 4000 1
## reject_deptB:applicant_gendermale 4000 1
## reject_deptC:applicant_gendermale 4000 1
## reject_deptD:applicant_gendermale 4000 1
## reject_deptE:applicant_gendermale 4000 1
## reject_deptF:applicant_gendermale 4000 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma(admit) 203.69 41.01 142.10 303.57 4000 1
## sigma(reject) 256.18 50.63 179.50 372.64 4000 1
## rescor(admit,reject) 0.61 0.17 0.21 0.85 4000 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(mQ1.brms3, ask = FALSE)
#(2)
mQ1.brms4 <- brm(admit_p ~ 0 + dept * applicant_gender,#separate intercept for each dept
prior=set_prior("normal(0,10)", class="b"),
data = d)
## Compiling the C++ model
## Start sampling
## Warning: There were 402 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low
## Warning: Examine the pairs() plot to diagnose sampling problems
summary(mQ1.brms4)
## Warning: The model has not converged (some Rhats are > 1.1). Do not analyse the results!
## We recommend running more iterations and/or setting stronger priors.
## Family: gaussian (identity)
## Formula: admit_p ~ 0 + dept * applicant_gender
## Data: d (Number of observations: 12)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## deptA 0.72 1.93 -3.75 4.79 2615
## deptB 0.69 2.14 -3.78 5.64 2925
## deptC 0.28 2.14 -4.33 4.74 3531
## deptD 0.40 2.25 -4.45 5.44 3637
## deptE 0.24 2.13 -4.45 4.93 3235
## deptF 0.13 2.16 -4.59 4.54 4000
## applicant_gendermale -0.08 2.19 -5.00 4.74 1754
## deptB:applicant_gendermale -0.03 3.38 -7.29 7.24 2246
## deptC:applicant_gendermale 0.16 3.37 -7.10 7.34 2234
## deptD:applicant_gendermale -0.04 3.46 -7.48 7.75 2322
## deptE:applicant_gendermale 0.12 3.29 -6.89 7.55 2200
## deptF:applicant_gendermale 0.05 3.60 -7.62 8.25 2106
## Rhat
## deptA 1
## deptB 1
## deptC 1
## deptD 1
## deptE 1
## deptF 1
## applicant_gendermale 1
## deptB:applicant_gendermale 1
## deptC:applicant_gendermale 1
## deptD:applicant_gendermale 1
## deptE:applicant_gendermale 1
## deptF:applicant_gendermale 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 1.89 1.57 0.21 5.71 22 1.16
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(mQ1.brms4, ask = FALSE)
#Test whether the coefficient for gender is different from 0 in the brms model
hypothesis(mQ1.brms3,"admit_applicant_gendermale = 0")
## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (admit_applicant_... = 0 2.19 9.69 -16.89 20.69 NA
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.
plot(hypothesis(mQ1.brms3,"admit_applicant_gendermale = 0"))
hypothesis(mQ1.brms4,"applicant_gendermale = 0")
## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (applicant_gender... = 0 -0.08 2.19 -5 4.74 NA
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.
plot(hypothesis(mQ1.brms4,"applicant_gendermale = 0"))
###Q2:Refit models 12.1 and 12.2 (Rcode 12.2 and 12.3) with brms
data(reedfrogs)
d2 <- reedfrogs
d2
## density pred size surv propsurv
## 1 10 no big 9 0.9000000
## 2 10 no big 10 1.0000000
## 3 10 no big 7 0.7000000
## 4 10 no big 10 1.0000000
## 5 10 no small 9 0.9000000
## 6 10 no small 9 0.9000000
## 7 10 no small 10 1.0000000
## 8 10 no small 9 0.9000000
## 9 10 pred big 4 0.4000000
## 10 10 pred big 9 0.9000000
## 11 10 pred big 7 0.7000000
## 12 10 pred big 6 0.6000000
## 13 10 pred small 7 0.7000000
## 14 10 pred small 5 0.5000000
## 15 10 pred small 9 0.9000000
## 16 10 pred small 9 0.9000000
## 17 25 no big 24 0.9600000
## 18 25 no big 23 0.9200000
## 19 25 no big 22 0.8800000
## 20 25 no big 25 1.0000000
## 21 25 no small 23 0.9200000
## 22 25 no small 23 0.9200000
## 23 25 no small 23 0.9200000
## 24 25 no small 21 0.8400000
## 25 25 pred big 6 0.2400000
## 26 25 pred big 13 0.5200000
## 27 25 pred big 4 0.1600000
## 28 25 pred big 9 0.3600000
## 29 25 pred small 13 0.5200000
## 30 25 pred small 20 0.8000000
## 31 25 pred small 8 0.3200000
## 32 25 pred small 10 0.4000000
## 33 35 no big 34 0.9714286
## 34 35 no big 33 0.9428571
## 35 35 no big 33 0.9428571
## 36 35 no big 31 0.8857143
## 37 35 no small 31 0.8857143
## 38 35 no small 35 1.0000000
## 39 35 no small 33 0.9428571
## 40 35 no small 32 0.9142857
## 41 35 pred big 4 0.1142857
## 42 35 pred big 12 0.3428571
## 43 35 pred big 13 0.3714286
## 44 35 pred big 14 0.4000000
## 45 35 pred small 22 0.6285714
## 46 35 pred small 12 0.3428571
## 47 35 pred small 31 0.8857143
## 48 35 pred small 17 0.4857143
# make the tank cluster variable
d2$tank <- 1:nrow(d2)
#fit the model with stan
m12.1.stan <- map2stan(
alist(
surv ~ dbinom( density , p ) ,
logit(p) <- a_tank[tank] ,
a_tank[tank] ~ dnorm( 0 , 5 )
), data=d2 )
## In file included from file10a46fa2057e.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/var.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/tools/config.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config.hpp:39:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config/compiler/clang.hpp:196:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command line>:6:9: note: previous definition is here
## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1
## ^
## 1 warning generated.
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.174301 seconds (Warm-up)
## 0.154378 seconds (Sampling)
## 0.328679 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 0.0001 seconds (Sampling)
## 0.000103 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
precis(m12.1.stan,depth = 2)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_tank[1] 2.54 1.26 0.73 4.36 484 1.00
## a_tank[2] 5.61 2.64 1.62 9.53 586 1.00
## a_tank[3] 0.89 0.72 -0.19 2.00 771 1.00
## a_tank[4] 5.73 2.79 1.63 9.75 666 1.00
## a_tank[5] 2.54 1.23 0.74 4.42 764 1.00
## a_tank[6] 2.58 1.28 0.76 4.50 639 1.00
## a_tank[7] 5.78 2.87 1.85 9.83 407 1.01
## a_tank[8] 2.59 1.16 0.78 4.27 728 1.00
## a_tank[9] -0.45 0.70 -1.61 0.62 1000 1.00
## a_tank[10] 2.57 1.21 0.69 4.37 572 1.00
## a_tank[11] 0.94 0.74 -0.31 2.01 818 1.00
## a_tank[12] 0.44 0.69 -0.60 1.54 1000 1.00
## a_tank[13] 0.93 0.75 -0.34 1.95 1000 1.00
## a_tank[14] 0.01 0.68 -1.08 1.05 1000 1.00
## a_tank[15] 2.48 1.10 0.75 4.06 919 1.00
## a_tank[16] 2.52 1.18 0.80 4.37 481 1.00
## a_tank[17] 3.53 1.21 1.70 5.29 559 1.00
## a_tank[18] 2.56 0.75 1.46 3.73 1000 1.00
## a_tank[19] 2.09 0.65 1.07 3.06 1000 1.00
## a_tank[20] 6.54 2.81 2.50 10.28 510 1.00
## a_tank[21] 2.62 0.81 1.48 3.96 887 1.00
## a_tank[22] 2.64 0.77 1.42 3.78 720 1.00
## a_tank[23] 2.59 0.77 1.51 3.91 785 1.00
## a_tank[24] 1.75 0.59 0.87 2.64 814 1.00
## a_tank[25] -1.19 0.47 -1.91 -0.44 974 1.00
## a_tank[26] 0.10 0.40 -0.52 0.75 1000 1.00
## a_tank[27] -1.70 0.52 -2.51 -0.86 1000 1.00
## a_tank[28] -0.60 0.41 -1.25 0.03 1000 1.00
## a_tank[29] 0.07 0.42 -0.56 0.74 645 1.00
## a_tank[30] 1.41 0.51 0.62 2.18 1000 1.00
## a_tank[31] -0.78 0.45 -1.53 -0.11 920 1.00
## a_tank[32] -0.43 0.41 -1.02 0.32 1000 1.00
## a_tank[33] 3.77 1.08 2.00 5.21 604 1.00
## a_tank[34] 2.94 0.76 1.78 4.08 875 1.00
## a_tank[35] 2.98 0.77 1.80 4.16 866 1.00
## a_tank[36] 2.14 0.56 1.26 2.97 794 1.00
## a_tank[37] 2.13 0.54 1.29 2.96 865 1.00
## a_tank[38] 6.78 2.48 3.13 10.32 652 1.00
## a_tank[39] 2.95 0.80 1.60 4.05 986 1.00
## a_tank[40] 2.48 0.59 1.47 3.31 1000 1.00
## a_tank[41] -2.11 0.52 -2.97 -1.33 808 1.00
## a_tank[42] -0.67 0.36 -1.21 -0.08 1000 1.00
## a_tank[43] -0.54 0.36 -1.07 0.10 1000 1.00
## a_tank[44] -0.43 0.34 -0.97 0.12 885 1.00
## a_tank[45] 0.54 0.35 -0.03 1.08 1000 1.00
## a_tank[46] -0.68 0.36 -1.22 -0.10 672 1.00
## a_tank[47] 2.16 0.57 1.26 3.06 794 1.00
## a_tank[48] -0.06 0.32 -0.59 0.43 1000 1.00
#refit the model 12.1 with brms
#(1)without interaction
m12.1.brms1 <- brm(propsurv ~ pred + size,
prior=set_prior("normal(0,5)", class="b"),
data = d2)
## Compiling the C++ model
## Start sampling
summary(m12.1.brms1)
## Family: gaussian (identity)
## Formula: propsurv ~ pred + size
## Data: d2 (Number of observations: 48)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 0.88 0.04 0.79 0.96 4000 1
## predpred -0.40 0.05 -0.50 -0.31 4000 1
## sizesmall 0.09 0.05 -0.01 0.19 4000 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 0.17 0.02 0.14 0.22 4000 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(m12.1.brms1,ask=FALSE)
#(2)with
m12.1.brms2 <- brm(propsurv ~ pred*size,
prior=set_prior("normal(0,5)", class="b"),
data = d2)
## Compiling the C++ model
## Start sampling
summary(m12.1.brms2)
## Family: gaussian (identity)
## Formula: propsurv ~ pred * size
## Data: d2 (Number of observations: 48)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 0.93 0.05 0.83 1.02 2793 1
## predpred -0.50 0.07 -0.64 -0.37 2522 1
## sizesmall -0.01 0.07 -0.15 0.13 2449 1
## predpred:sizesmall 0.20 0.10 0.01 0.40 2250 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 0.17 0.02 0.14 0.21 3671 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(m12.1.brms2,ask=FALSE)
#fit model 12.2 with stan
m12.2.stan <- map2stan(
alist(
surv ~ dbinom( density , p ) ,
logit(p) <- a_tank[tank] ,
a_tank[tank] ~ dnorm( a , sigma ) ,
a ~ dnorm(0,1) ,
sigma ~ dcauchy(0,1)
), data=d2 , chains=4 )
## In file included from file10a478873559.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/var.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/tools/config.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config.hpp:39:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config/compiler/clang.hpp:196:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command line>:6:9: note: previous definition is here
## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1
## ^
## 1 warning generated.
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.177463 seconds (Warm-up)
## 0.105072 seconds (Sampling)
## 0.282535 seconds (Total)
## The following numerical problems occured the indicated number of times on chain 1
## count
## Exception thrown at line 17: normal_log: Scale parameter is 0, but must be > 0! 1
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.18435 seconds (Warm-up)
## 0.112896 seconds (Sampling)
## 0.297246 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.174656 seconds (Warm-up)
## 0.103942 seconds (Sampling)
## 0.278598 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 4).
##
## Chain 4, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.180502 seconds (Warm-up)
## 0.107807 seconds (Sampling)
## 0.288309 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 8.1e-05 seconds (Sampling)
## 8.5e-05 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
precis(m12.2.stan, depth = 2)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_tank[1] 2.12 0.87 0.84 3.60 4000 1
## a_tank[2] 3.04 1.10 1.33 4.72 4000 1
## a_tank[3] 1.00 0.67 -0.07 2.02 4000 1
## a_tank[4] 3.04 1.11 1.21 4.66 4000 1
## a_tank[5] 2.13 0.89 0.59 3.44 4000 1
## a_tank[6] 2.11 0.86 0.74 3.40 4000 1
## a_tank[7] 3.03 1.10 1.34 4.70 4000 1
## a_tank[8] 2.12 0.84 0.83 3.46 4000 1
## a_tank[9] -0.19 0.63 -1.20 0.80 4000 1
## a_tank[10] 2.10 0.85 0.72 3.44 4000 1
## a_tank[11] 1.00 0.67 -0.05 2.06 4000 1
## a_tank[12] 0.57 0.61 -0.35 1.58 4000 1
## a_tank[13] 1.00 0.67 0.00 2.13 4000 1
## a_tank[14] 0.19 0.62 -0.78 1.19 4000 1
## a_tank[15] 2.13 0.87 0.78 3.47 4000 1
## a_tank[16] 2.11 0.84 0.74 3.40 4000 1
## a_tank[17] 2.89 0.78 1.68 4.10 4000 1
## a_tank[18] 2.40 0.68 1.25 3.37 4000 1
## a_tank[19] 2.00 0.57 1.05 2.83 4000 1
## a_tank[20] 3.67 1.00 2.15 5.16 4000 1
## a_tank[21] 2.39 0.67 1.35 3.45 4000 1
## a_tank[22] 2.38 0.66 1.40 3.47 4000 1
## a_tank[23] 2.39 0.66 1.34 3.42 4000 1
## a_tank[24] 1.70 0.54 0.79 2.49 4000 1
## a_tank[25] -1.00 0.45 -1.69 -0.27 4000 1
## a_tank[26] 0.15 0.40 -0.47 0.82 4000 1
## a_tank[27] -1.43 0.49 -2.19 -0.66 4000 1
## a_tank[28] -0.47 0.41 -1.08 0.21 4000 1
## a_tank[29] 0.15 0.41 -0.51 0.80 4000 1
## a_tank[30] 1.44 0.49 0.73 2.28 4000 1
## a_tank[31] -0.64 0.42 -1.29 0.06 4000 1
## a_tank[32] -0.31 0.41 -0.91 0.37 4000 1
## a_tank[33] 3.18 0.78 2.01 4.43 4000 1
## a_tank[34] 2.71 0.65 1.73 3.75 4000 1
## a_tank[35] 2.70 0.66 1.67 3.72 4000 1
## a_tank[36] 2.06 0.51 1.23 2.82 4000 1
## a_tank[37] 2.05 0.49 1.28 2.83 4000 1
## a_tank[38] 3.86 0.96 2.29 5.24 4000 1
## a_tank[39] 2.71 0.67 1.65 3.76 4000 1
## a_tank[40] 2.34 0.56 1.43 3.16 4000 1
## a_tank[41] -1.82 0.48 -2.61 -1.11 4000 1
## a_tank[42] -0.56 0.35 -1.16 -0.03 4000 1
## a_tank[43] -0.46 0.34 -0.98 0.08 4000 1
## a_tank[44] -0.33 0.34 -0.86 0.23 4000 1
## a_tank[45] 0.59 0.34 0.04 1.12 4000 1
## a_tank[46] -0.57 0.35 -1.10 0.00 4000 1
## a_tank[47] 2.06 0.52 1.22 2.84 4000 1
## a_tank[48] 0.00 0.34 -0.55 0.52 4000 1
## a 1.30 0.26 0.88 1.69 4000 1
## sigma 1.62 0.22 1.29 1.96 2729 1
#refit the model 12.2 with brms, by specifying different priors for specific coefficients
#(1)without interaction
m12.2.brms1 <- brm(propsurv ~ pred + size,
data = d2,
prior= c(
set_prior("normal(0,1)", class="Intercept"),
set_prior("normal(0,1)", class="b"),
set_prior("cauchy(0,1)", class = "sigma")))
## Compiling the C++ model
## Start sampling
summary(m12.2.brms1)
## Family: gaussian (identity)
## Formula: propsurv ~ pred + size
## Data: d2 (Number of observations: 48)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 0.88 0.04 0.79 0.96 4000 1
## predpred -0.40 0.05 -0.50 -0.30 4000 1
## sizesmall 0.09 0.05 -0.01 0.19 4000 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 0.17 0.02 0.14 0.22 4000 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(m12.2.brms1,ask=FALSE)
#(2)with interaction
m12.2.brms2 <- brm(propsurv ~ pred * size,
data = d2,
prior= c(
set_prior("normal(0,1)", class="Intercept"),
set_prior("normal(0,1)", class="b"),
set_prior("cauchy(0,1)", class = "sigma")))
## Compiling the C++ model
## Start sampling
summary(m12.2.brms2)
## Family: gaussian (identity)
## Formula: propsurv ~ pred * size
## Data: d2 (Number of observations: 48)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 0.92 0.05 0.83 1.02 2459 1
## predpred -0.50 0.07 -0.63 -0.37 2435 1
## sizesmall 0.00 0.07 -0.14 0.13 2264 1
## predpred:sizesmall 0.19 0.10 0.00 0.38 2175 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 0.17 0.02 0.14 0.21 3294 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(m12.2.brms2,ask=FALSE)
#comparing models
loo(m12.1.brms1,m12.1.brms2,m12.2.brms1,m12.2.brms2)
## Warning: Some Pareto k diagnostic values are slightly high. See
## help('pareto-k-diagnostic') for details.
## Warning: Some Pareto k diagnostic values are slightly high. See
## help('pareto-k-diagnostic') for details.
## LOOIC SE
## m12.1.brms1 -28.91 10.61
## m12.1.brms2 -30.49 12.27
## m12.2.brms1 -28.84 10.68
## m12.2.brms2 -30.48 12.24
## m12.1.brms1 - m12.1.brms2 1.58 4.69
## m12.1.brms1 - m12.2.brms1 -0.07 0.10
## m12.1.brms1 - m12.2.brms2 1.57 4.50
## m12.1.brms2 - m12.2.brms1 -1.65 4.74
## m12.1.brms2 - m12.2.brms2 -0.01 0.23
## m12.2.brms1 - m12.2.brms2 1.65 4.55
d3 <- read.csv("/Users/xyyue/Documents/RClub/Rclub-rethinking_Xiaoyan.Yue/Assignment_Chapter_09/TomatoR2CSHL.csv")
head(d3)
## shelf flat col row acs trt days date hyp int1 int2 int3 int4
## 1 Z 1 B 1 LA2580 H 28 5/5/08 19.46 2.37 1.59 1.87 0.51
## 2 Z 1 C 1 LA1305 H 28 5/5/08 31.28 3.34 0.01 9.19 1.62
## 3 Z 1 D 1 LA1973 H 28 5/5/08 56.65 8.43 2.39 6.70 3.69
## 4 Z 1 E 1 LA2748 H 28 5/5/08 35.18 0.56 0.00 1.60 0.61
## 5 Z 1 F 1 LA2931 H 28 5/5/08 35.32 0.82 0.02 1.49 0.46
## 6 Z 1 G 1 LA1317 H 28 5/5/08 28.74 1.07 6.69 5.72 4.76
## intleng totleng petleng leafleng leafwid leafnum ndvi lat lon
## 1 6.34 25.80 15.78 30.53 34.44 5 111 -9.5167 -78.0083
## 2 14.16 45.44 12.36 22.93 13.99 4 120 -13.3833 -75.3583
## 3 21.21 77.86 13.05 46.71 43.78 5 110 -16.2333 -71.7000
## 4 2.77 37.95 8.08 26.82 33.28 5 105 -20.4833 -69.9833
## 5 2.79 38.11 7.68 22.40 23.61 5 106 -20.9167 -69.0667
## 6 18.24 46.98 23.66 42.35 42.35 5 132 -13.4167 -73.8417
## alt species who
## 1 740 S. pennellii Dan
## 2 3360 S. peruvianum Dan
## 3 2585 S. peruvianum Dan
## 4 1020 S. chilense Dan
## 5 2460 S. chilense Dan
## 6 2000 S. chmielewskii Dan
#fit a model with stan for intleng as a function of species, trt and their interaction
#categorical variable for trt, species and shelf
d3$trt2 <- as.numeric(d3$trt)-1 #0 = H, 1=L
d3$species_id <- coerce_index(d3$species)
d3$shelf_id <- coerce_index(d3$shelf)
d3.subset <- d3[,c("shelf_id","trt2","intleng","species_id")]
head(d3.subset)
## shelf_id trt2 intleng species_id
## 1 6 0 6.34 4
## 2 6 0 14.16 5
## 3 6 0 21.21 5
## 4 6 0 2.77 1
## 5 6 0 2.79 1
## 6 6 0 18.24 2
summary(d3.subset)
## shelf_id trt2 intleng species_id
## Min. :1.000 Min. :0.0000 Min. : 0.000 Min. :1.000
## 1st Qu.:2.000 1st Qu.:0.0000 1st Qu.: 9.637 1st Qu.:2.000
## Median :3.000 Median :1.0000 Median :17.255 Median :3.000
## Mean :3.512 Mean :0.5089 Mean :20.340 Mean :2.927
## 3rd Qu.:5.000 3rd Qu.:1.0000 3rd Qu.:28.145 3rd Qu.:4.000
## Max. :6.000 Max. :1.0000 Max. :92.420 Max. :5.000
#fit the model
mQ3.stan <- map2stan(
alist(
intleng ~ dnorm(mu,sigma),
mu <- a + a_shelf[shelf_id] + bs[species_id] + bt*trt2 + bst[species_id]*trt2,
a ~ dnorm(20,100),
a_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
bs[species_id] ~ dnorm(0,10),
bt ~ dnorm(0,20),
bst[species_id] ~ dnorm(0,10),
sigma ~ dcauchy(0,1),
sigma_shelf ~ dexp(1)
),
data = d3.subset,
iter = 4000,
warmup = 1000,
chains=4,
cores=2
)
## In file included from file10a42593103f.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/var.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/tools/config.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config.hpp:39:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config/compiler/clang.hpp:196:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command line>:6:9: note: previous definition is here
## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1
## ^
## 1 warning generated.
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
##
## SAMPLING FOR MODEL 'intleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 0.000257 seconds (Sampling)
## 0.000261 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]
## Warning in map2stan(alist(intleng ~ dnorm(mu, sigma), mu <- a + a_shelf[shelf_id] + : There were 1 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
precis(mQ3.stan)
## Warning in precis(mQ3.stan): There were 1 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
## 16 vector or matrix parameters omitted in display. Use depth=2 to show them.
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 10.67 4.62 2.93 17.67 3228 1
## bt 17.03 4.64 9.27 24.07 3461 1
## sigma 10.61 0.24 10.20 10.97 10429 1
## sigma_shelf 1.42 0.68 0.39 2.44 2920 1
pairs(mQ3.stan)
#fit a model with brm
#(1)without interaction
mQ3.brms1 <- brm(intleng ~ 0 + species + trt + (1|shelf),
data = d3,
prior = c(
set_prior("normal(0,10)",class="b"), # sets prior for all b coefficients
set_prior("normal(0,5)",class="b", coef = "trtL"), #set prior for "trtL"
set_prior("cauchy(0,1)", class = "sigma"), #half cauchy prior for sigma
set_prior("normal(0,1)", class = "sd", group = "shelf") #prior for shelf
)
)
## Compiling the C++ model
## Start sampling
summary(mQ3.brms1)
## Family: gaussian (identity)
## Formula: intleng ~ 0 + species + trt + (1 | shelf)
## Data: d3 (Number of observations: 1008)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Group-Level Effects:
## ~shelf (Number of levels: 6)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 1.29 0.49 0.44 2.36 1460 1
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## speciesS.chilense 7.82 1.11 5.57 9.98 1561 1
## speciesS.chmielewskii 11.89 1.10 9.68 14.07 1566 1
## speciesS.habrochaites 15.11 1.11 13.00 17.38 1579 1
## speciesS.pennellii 5.49 1.24 3.13 7.93 1921 1
## speciesS.peruvianum 13.75 1.11 11.54 15.94 1487 1
## trtL 17.13 1.30 14.30 19.50 1302 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 10.67 0.24 10.21 11.13 4000 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(mQ3.brms1)
#(2)with interaction
mQ3.brms2 <- brm(intleng ~ 0 + trt * species + (1|shelf),
data = d3,
prior = c(
set_prior("normal(0,10)",class="b"), # sets prior for all b coefficients
set_prior("normal(0,5)",class="b", coef = "trtL"), #set prior for "trtL"
set_prior("cauchy(0,1)", class = "sigma"), #half cauchy prior for sigma
set_prior("normal(0,1)", class = "sd", group = "shelf") #prior for shelf
)
)
## Compiling the C++ model
## Start sampling
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
summary(mQ3.brms2)
## Family: gaussian (identity)
## Formula: intleng ~ 0 + trt * species + (1 | shelf)
## Data: d3 (Number of observations: 1008)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Group-Level Effects:
## ~shelf (Number of levels: 6)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 1.39 0.56 0.39 2.62 1270 1
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## trtH 5.87 1.30 3.24 8.39 1745
## trtL 25.57 1.41 22.61 28.13 1077
## speciesS.chmielewskii 7.75 1.41 4.93 10.45 2120
## speciesS.habrochaites 9.30 1.37 6.59 12.09 1925
## speciesS.pennellii -0.48 1.67 -3.80 2.74 2110
## speciesS.peruvianum 7.41 1.38 4.67 10.01 1817
## trtL:speciesS.chmielewskii -6.26 1.94 -10.05 -2.47 2141
## trtL:speciesS.habrochaites -3.04 1.92 -6.81 0.64 2114
## trtL:speciesS.pennellii -2.63 2.25 -7.09 1.73 1869
## trtL:speciesS.peruvianum -1.97 1.97 -5.85 1.91 2018
## Rhat
## trtH 1
## trtL 1
## speciesS.chmielewskii 1
## speciesS.habrochaites 1
## speciesS.pennellii 1
## speciesS.peruvianum 1
## trtL:speciesS.chmielewskii 1
## trtL:speciesS.habrochaites 1
## trtL:speciesS.pennellii 1
## trtL:speciesS.peruvianum 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 10.61 0.24 10.16 11.07 4000 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(mQ3.brms2)
#compare the model
loo(mQ3.brms1, mQ3.brms2)
## LOOIC SE
## mQ3.brms1 7645.14 64.73
## mQ3.brms2 7638.18 65.98
## mQ3.brms1 - mQ3.brms2 6.96 5.60